/*
 * stack_description.c
 *
 *  Created on: Mar 12, 2012
 *      Author: qian
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "stack.h"
#include "auxiliary.h"


/*Construct stack description from layer configuration file and floorplan*/
void build_stack_description(IO_t *IO, stk_descr_t *stk)
{
	FILE *fp=fopen(IO->lcf, "r");
	parse_config_file(stk, fp);
	set_bgmap(stk);

}

/* parse the layer file open for reading	*/
void parse_config_file(stk_descr_t *stk, FILE *fp)
{
	char line[LINE_SIZE], *ptr, cval;
	char line_cpy[LINE_SIZE];
	int count, i = 0, field = LCF_DIGIT, ival, ival1;
	double dval, dval1, dval2, dval3;

	fseek(fp, 0, SEEK_SET);

	/*read the line for total layer number and cooling method*/
	count = 0;

	while(count<1)
	{
		fgets(line, LINE_SIZE, fp);
		if (feof(fp))
		{
			fatal("Please Check CONFIG file format.");
		}
		memcpy(line_cpy, line, LINE_SIZE);
		/* ignore comments and empty lines	*/
		ptr = strtok(line_cpy, " \r\t\n");
		if (!ptr || ptr[0] == '#')
			continue;
		if(sscanf(line, "%d%d", &ival, &ival1)!=2)
			fatal("invalid total layer number and cooling method.");
		stk->n_layers=ival;
		stk->cooling_method=ival1;

		count++;
	}
	stk->layers=(layer_t*)calloc(stk->n_layers,sizeof(layer_t));

	/*read every circuit layer's property*/
	count = 0;
	while (count < (stk->n_layers * LCF_NPARAMS)) {
		fgets(line, LINE_SIZE, fp);
		if (feof(fp))
		{
			fatal("Please Check layer spec section of CONFIG file format.");
		}

		memcpy(line_cpy, line, LINE_SIZE);
		/* ignore comments and empty lines	*/
		ptr = strtok(line_cpy, " \r\t\n");
		if (!ptr || ptr[0] == '#')
			continue;

		switch (field)
		{
			case LCF_DIGIT:
						if (sscanf(line, "%d%lf%lf%lf", &ival, &dval, &dval1, &dval2) != 4)
							fatal("invalid DIGIT line describing layer number, thickness of BEOL, SOURCE and BULK\n");
						if(ival >= stk->n_layers || ival < 0)
							fatal("layer number must be >= 0 and < no. of layers\n");
						i = ival;
						stk->layers[i].idx = ival;
						stk->layers[i].BEOL_height=dval;
						stk->layers[i].source_height=dval1;
						stk->layers[i].bulk_height=dval2;
						break;

			case LCF_STRING:
						stk->layers[i].flp = read_flp(ptr);
						/* first layer	*/
						if (count < LCF_NPARAMS)
						{
							stk->width = get_total_width(stk->layers[i].flp);
							stk->length = get_total_length(stk->layers[i].flp);
						}
						else if(!eq(stk->width, get_total_width(stk->layers[i].flp)) ||
								!eq(stk->length, get_total_length(stk->layers[i].flp)))
							fatal("width and height differ across layers\n");
						break;
			default:
						fatal("invalid field id\n");
						break;
		}
		field = (field + 1) % LCF_NPARAMS;
		count++;
	}


	/*read cooling package information*/
	count=0;
	if(stk->cooling_method==1) //heatsink cooling
	{
		stk->heatsink=(heatsink_t *)malloc(sizeof(heatsink_t));

		while(count<1) //heatsink specs
		{
			fgets(line, LINE_SIZE, fp);
			if (feof(fp))
			{
				printf("unexpected end of file! Please check heatsink section of CONFIG file");
				break;
			}
			memcpy(line_cpy, line, LINE_SIZE);
			/* ignore comments and empty lines	*/
			ptr = strtok(line_cpy, " \r\t\n");
			if (!ptr || ptr[0] == '#')
				continue;

			if(sscanf(line, "%lf%lf%lf%lf", &dval, &dval1, &dval2, &dval3)!=4)
				printf("invalid line for heatsink spec");

			stk->heatsink->t_heatsink=dval;
			stk->heatsink->width_heatsink=dval1;
			stk->heatsink->k_heatsink=dval2;
			stk->heatsink->r_conv=dval3;

			count++;
		}

		while(count<2) //spreader specs
		{
			fgets(line, LINE_SIZE, fp);
			if (feof(fp))
			{
				printf("unexpected end of file! Please check spreader section of CONFIG file");
				break;
			}
			memcpy(line_cpy, line, LINE_SIZE);
			/* ignore comments and empty lines	*/
			ptr = strtok(line_cpy, " \r\t\n");
			if (!ptr || ptr[0] == '#')
				continue;

			if(sscanf(line, "%lf%lf%lf", &dval, &dval1, &dval2)!=3)
				printf("invalid line for spreader spec");

			stk->heatsink->t_spreader=dval;
			stk->heatsink->width_spreader=dval1;
			stk->heatsink->k_spreader=dval2;

			count++;
		}


		while(count<3) //TIM specs
		{
			fgets(line, LINE_SIZE, fp);
			if (feof(fp))
			{
				printf("unexpected end of file! Please check TIM section of CONFIG file");
				break;
			}
			memcpy(line_cpy, line, LINE_SIZE);
			/* ignore comments and empty lines	*/
			ptr = strtok(line_cpy, " \r\t\n");
			if (!ptr || ptr[0] == '#')
				continue;

			if(sscanf(line, "%lf%lf", &dval, &dval1)!=2)
				printf("invalid line for TIM spec");

			stk->heatsink->t_TIM=dval;
			stk->heatsink->k_TIM=dval1;

			count++;
		}

	}
	else if(stk->cooling_method==2) //microfluidic cooling
	{
		stk->microchannel=(microchannel_t *)malloc(sizeof(microchannel_t));

		while(count<1)
		{
			fgets(line, LINE_SIZE, fp);
			if (feof(fp))
			{
				printf("unexpected end of file! Please check microchannel section of CONFIG file");
				break;
			}
			memcpy(line_cpy, line, LINE_SIZE);
			/* ignore comments and empty lines	*/
			ptr = strtok(line_cpy, " \r\t\n");
			if (!ptr || ptr[0] == '#')
				continue;

			if(sscanf(line, "%lf%lf%lf%lf", &dval, &dval1, &dval2, &dval3)!=4)
				fatal("invalid line for microchannel spec");

			stk->microchannel->channel_width=dval;
			stk->microchannel->channel_height=dval1;
			stk->microchannel->wall_width=dval2;
			stk->microchannel->default_flowrate=dval3;

			count++;
		}

		stk->microchannel->n_channels_per_layer=double2int(stk->width/(stk->microchannel->channel_width+stk->microchannel->wall_width));
	}

	/*read grid resolution*/
	count=0;
	while(count<1)
	{
		fgets(line, LINE_SIZE, fp);
		if (feof(fp))
		{
			printf("unexpected end of file! Please check grid resolution section of CONFIG file");
			break;
		}
		memcpy(line_cpy, line, LINE_SIZE);
		/* ignore comments and empty lines	*/
		ptr = strtok(line_cpy, " \r\t\n");
		if (!ptr || ptr[0] == '#')
			continue;

		if(sscanf(line, "%d%d", &ival, &ival1)!=2)
			fatal("invalid line for grid resolution");

		stk->rows=ival;
		stk->cols=ival1;

		count++;
	}

	/*TSV map section*/
	count=0;

	while(count<1)
	{
		fgets(line, LINE_SIZE, fp);
		if (feof(fp))
		{
			printf("unexpected end of file! Please check TSV section of CONFIG file");
			break;
		}
		ptr = strtok(line, " \r\t\n");
		if (!ptr || ptr[0] == '#')
			continue;

		if(strcmp(ptr,"null")!=0)
		{
			read_tsv(stk, ptr);
			tsv_translate(stk);
		}
		else
		{
			printf("TSV file not specified, ignoring effect of TSVs\n");
			init_tsv(stk);
			stk->tsv_physi=(tsv_physi_t *)malloc(sizeof(tsv_physi_t));
		}

		count++;
	}


	/* allocate the block-grid maps */
		for(i=0; i < stk->n_layers; i++) {
			stk->layers[i].b2gmap = new_b2gmap(stk->rows, stk->cols);
			stk->layers[i].g2bmap = (glist_t *) calloc(stk->layers[i].flp->n_units,
									   sizeof(glist_t));
			if (!stk->layers[i].g2bmap)
				fatal("memory allocation error0\n");
		}

}

void read_tsv(stk_descr_t *stk, char *file)
{
	char str[STR_SIZE];
	FILE *fp;
	int ival, ival1, ival2, ival3, count;
	double dval, dval1, dval2;
	char *ptr, line[LINE_SIZE], line_cpy[LINE_SIZE];

	fp=fopen(file, "r");
	if (!fp) {
		sprintf(str, "error opening tsv config file %s\n", file);
		fatal(str);
	}

	fseek(fp, 0, SEEK_SET);

	/*read the line for total TSV number and material's thermal conductivity*/
	count = 0;

	while(count<1)
	{
		fgets(line, LINE_SIZE, fp);
		if (feof(fp))
		{
			fatal("Please Check TSV config file format.");
		}
		memcpy(line_cpy, line, LINE_SIZE);
		/* ignore comments and empty lines	*/
		ptr = strtok(line_cpy, " \r\t\n");
		if (!ptr || ptr[0] == '#')
			continue;
		if(sscanf(line, "%d%lf", &ival, &dval)!=2)
			fatal("invalid line for total TSV number");

		stk->n_tsvs=ival;
		stk->k_tsv=dval;
		stk->tsv_physi=(tsv_physi_t *)malloc(ival*sizeof(tsv_physi_t));

		count++;
	}

	/*read TSV positions and sizes*/

	count=0;

	while(count<stk->n_tsvs)
	{
		fgets(line, LINE_SIZE, fp);
		if (feof(fp))
		{
			fatal("Please Check TSV config file format.");
		}
		memcpy(line_cpy, line, LINE_SIZE);
		/* ignore comments and empty lines	*/
		ptr = strtok(line_cpy, " \r\t\n");
		if (!ptr || ptr[0] == '#')
			continue;

		if(sscanf(line, "%lf%lf%d%d%d%d%lf", &dval, &dval1, &ival, &ival1, &ival2, &ival3, &dval2)!=7)
			fatal("invalid line for TSV properties");

		stk->tsv_physi[count].x=dval;
		stk->tsv_physi[count].y=dval1;
		stk->tsv_physi[count].z1=ival;
		stk->tsv_physi[count].sl1=ival1;
		stk->tsv_physi[count].z2=ival2;
		stk->tsv_physi[count].sl2=ival3;
		stk->tsv_physi[count].width=dval2;

		count++;
	}


}

/*parse floorplan files*/
flp_t *read_flp(char *file)
{
	char str[STR_SIZE];
	FILE *fp;
	flp_t *flp;
	int count, i, j;

	fp = fopen (file, "r");

	if (!fp) {
		sprintf(str, "error opening file %s\n", file);
		fatal(str);
	}

	/* 1st pass - find n_units	*/
	count = flp_count_units(fp);
	if(!count)
		fatal("no units specified in the floorplan file\n");

	/* allocate initial memory */
	flp = flp_alloc_init_mem(count);

	/* 2nd pass - populate block info	*/
	flp_populate_blks(flp, fp);

	fclose(fp);

	/* make sure the origin is (0,0)	*/
	flp_translate(flp, 0, 0);
	return flp;
}

/*
 * find the number of units from the
 * floorplan file
 */
int flp_count_units(FILE *fp)
{
    char str1[LINE_SIZE], str2[LINE_SIZE];
	char name[STR_SIZE];
	double leftx, bottomy, width, height, power;
	char *ptr;
    int count = 0;

	fseek(fp, 0, SEEK_SET);
	while(!feof(fp)) {
		fgets(str1, LINE_SIZE, fp);
		if (feof(fp))
			break;
		strcpy(str2, str1);

		/* ignore comments and empty lines	*/
		ptr = strtok(str1, " \r\t\n");
		if (!ptr || ptr[0] == '#')
			continue;

		/* functional block placement information	*/
		if (sscanf(str2, "%s%lf%lf%lf%lf%lf", name, &leftx, &bottomy,
		  		   &width, &height, &power) == 6)
			count++;
	}
	return count;
}

/*allocate memory for floorplan units*/
flp_t *flp_alloc_init_mem(int count)
{
	int i;
	flp_t *flp;
	flp = (flp_t *) calloc (1, sizeof(flp_t));
	if(!flp)
		fatal("memory allocation error1\n");
	flp->units = (unit_t *) calloc(count, sizeof(unit_t));
	if (!flp->units)
		fatal("memory allocation error2\n");
	flp->n_units = count;

	return flp;
}

/* populate block information	*/
void flp_populate_blks(flp_t *flp, FILE *fp)
{
	int i=0;
	char str[LINE_SIZE], copy[LINE_SIZE];
	char name1[STR_SIZE];
	double width, height, leftx, bottomy, power;
	char *ptr;

	fseek(fp, 0, SEEK_SET);
	while(!feof(fp)) {		/* second pass	*/
		fgets(str, LINE_SIZE, fp);
		if (feof(fp))
			break;
		strcpy(copy, str);

		/* ignore comments and empty lines	*/
		ptr = strtok(str, " \r\t\n");
		if (!ptr || ptr[0] == '#')
			continue;

		if (sscanf(copy, "%s%lf%lf%lf%lf%lf", name1, &width, &height,
				   &leftx, &bottomy, &power) == 6) {
			strcpy(flp->units[i].name, name1);
			flp->units[i].width = width;
			flp->units[i].height = height;
			flp->units[i].leftx = leftx;
			flp->units[i].bottomy = bottomy;
			flp->units[i].power=power;
			i++;
			/* skip connectivity info	*/
		}
	}
	if (i != flp->n_units)
	  fatal("mismatch of number of units\n");
}

/* translate the floorplan to new origin (x,y)	*/
void flp_translate(flp_t *flp, double x, double y)
{
	int i;
	double minx = flp->units[0].leftx;
	double miny = flp->units[0].bottomy;

	for (i=1; i < flp->n_units; i++) {
		if (minx > flp->units[i].leftx)
			minx = flp->units[i].leftx;
		if (miny > flp->units[i].bottomy)
			miny = flp->units[i].bottomy;
	}
	for (i=0; i < flp->n_units; i++) {
		flp->units[i].leftx += (x - minx);
		flp->units[i].bottomy += (y - miny);
	}
}

/*obtain die width from floorplan information*/
double get_total_width(flp_t *flp)
{
	int i;
	double min_x = flp->units[0].leftx;
	double max_x = flp->units[0].leftx + flp->units[0].width;

	for (i=1; i < flp->n_units; i++) {
		if (flp->units[i].leftx < min_x)
			min_x = flp->units[i].leftx;
		if (flp->units[i].leftx + flp->units[i].width > max_x)
			max_x = flp->units[i].leftx + flp->units[i].width;
	}

	return (max_x - min_x);
}

/*obtain die length from floorplan information*/
double get_total_length(flp_t *flp)
{
	int i;
	double min_y = flp->units[0].bottomy;
	double max_y = flp->units[0].bottomy + flp->units[0].height;

	for (i=1; i < flp->n_units; i++) {
		if (flp->units[i].bottomy < min_y)
			min_y = flp->units[i].bottomy;
		if (flp->units[i].bottomy + flp->units[i].height > max_y)
			max_y = flp->units[i].bottomy + flp->units[i].height;
	}

	return (max_y - min_y);
}

/*allocate memory for b2gmap*/
blist_t ***new_b2gmap(int rows, int cols)
{
	int i;
	blist_t ***b2gmap;

	b2gmap = (blist_t ***) calloc (rows, sizeof(blist_t **));
	b2gmap[0] = (blist_t **) calloc (rows * cols, sizeof(blist_t *));
	if (!b2gmap || !b2gmap[0])
		fatal("memory allocation error3\n");

	for(i=1; i < rows; i++)
		b2gmap[i] = b2gmap[0] + cols * i;

	return b2gmap;
}

/* destructor	*/
void delete_b2gmap(blist_t ***b2gmap, int rows, int cols)
{
	int i, j;
	blist_t *ptr, *temp;

	/* free the linked list	*/
	for(i=0; i < rows; i++)
		for(j=0; j < cols; j++) {
			ptr = b2gmap[i][j];
			while(ptr) {
				temp = ptr->next;
				free(ptr);
				ptr = temp;
			}
		}

	/* free the array space	*/
	free(b2gmap[0]);
	free(b2gmap);
}

/* compute the power/temperature average weighted by occupancies	*/
double blist_avg(blist_t *ptr, flp_t *flp, double *v, int type)
{
	double  val = 0.0;

	for(; ptr; ptr = ptr->next) {
		if (type == V_POWER)
		{
//			printf("idx is %d, width is %f, occupancy is %f, v is %f\n",
//					ptr->idx,flp->units[ptr->idx].width,ptr->occupancy, v[ptr->idx]);
			val += ptr->occupancy * v[ptr->idx] / (flp->units[ptr->idx].width *
				   flp->units[ptr->idx].height);
		}
		else if (type == V_TEMP)
			val += ptr->occupancy * v[ptr->idx];
		else
			fatal("unknown vector type\n");

	}
//	printf("val is %f\n", val);
	return val;
}

/* re-initialize */
void reset_b2gmap(stk_descr_t *stk)
{
	int n, i, j;
	blist_t *ptr, *temp;

	/* free the linked list	*/
	for(n=0;n<stk->n_layers;n++)
		for(i=0; i < stk->rows; i++)
			for(j=0; j < stk->cols; j++)
			{
				ptr = stk->layers[n].b2gmap[i][j];
				while(ptr) {
					temp = ptr->next;
					free(ptr);
					ptr = temp;
				}
				stk->layers[n].b2gmap[i][j] = NULL;
			}
}


/* create a linked list node and append it at the end	*/
void blist_append(blist_t *head, int idx, double occupancy)
{
	blist_t *tail = NULL;

	if(!head)
		fatal("blist_append called with empty list\n");

	/* traverse till the end	*/
	for(; head; head = head->next)
		tail = head;

	/* append */
	tail->next =  new_blist(idx, occupancy);
}


/* constructors	*/
blist_t *new_blist(int idx, double occupancy)
{
	blist_t *ptr = (blist_t *) calloc (1, sizeof(blist_t));
	if (!ptr)
		fatal("memory allocation error4\n");
	ptr->idx = idx;
	ptr->occupancy = occupancy;
	ptr->next = NULL;
	return ptr;
}

/* setup the block and grid mapping data structures	*/
void set_bgmap(stk_descr_t *stk)
{
	/* i1, i2, j1 and j2 are indices of the boundary grid cells	*/
	int n, i, j, u, i1, i2, j1, j2;
	layer_t *layer=stk->layers;

	/* shortcuts for cell width(cw) and cell height(ch)	*/
	double cw = stk->width / stk->cols;
	double ch = stk->length / stk->rows;

	/* initialize	*/
	reset_b2gmap(stk);

	/* for each functional unit	*/
	for(n=0;n<stk->n_layers;n++)
	for(u=0; u < layer[n].flp->n_units; u++) {
		/* shortcuts for unit boundaries	*/
		double lu = layer[n].flp->units[u].leftx;
		double ru = lu + layer[n].flp->units[u].width;
		double bu = layer[n].flp->units[u].bottomy;
		double tu = bu + layer[n].flp->units[u].height;

		/* top index (lesser row) = rows - ceil (topy / cell height)	*/
		i1 = stk->rows - tolerant_ceil(tu/ch);
		/* bottom index (greater row) = rows - floor (bottomy / cell height)	*/
		i2 = stk->rows - tolerant_floor(bu/ch);
		/* left index = floor (leftx / cell width)	*/
		j1 = tolerant_floor(lu/cw);
		/* right index = ceil (rightx / cell width)	*/
		j2 = tolerant_ceil(ru/cw);
		/* sanity check	*/
		if((i1 < 0) || (j1 < 0))
			fatal("negative grid cell start index!\n");
		if((i2 > stk->rows) || (j2 >stk->cols))
			fatal("grid cell end index out of bounds!\n");
		if((i1 >= i2) || (j1 >= j2))
			fatal("invalid floorplan spec or grid resolution\n");

		/* setup g2bmap	*/
		layer[n].g2bmap[u].i1 = i1;
		layer[n].g2bmap[u].i2 = i2;
		layer[n].g2bmap[u].j1 = j1;
		layer[n].g2bmap[u].j2 = j2;

		/* setup b2gmap	*/
		/* for each grid cell in this unit	*/
		for(i=i1; i < i2; i++)
			for(j=j1; j < j2; j++)
				/* grid cells fully overlapped by this unit	*/
				if ((i > i1) && (i < i2-1) && (j > j1) && (j < j2-1)) {
					/* first unit in the list	*/
					if (!layer[n].b2gmap[i][j])
						layer[n].b2gmap[i][j] = new_blist(u, 1.0);
					else {
					/* this should not occur since the grid cell is
					 * fully covered and hence, no other unit should
					 * be sharing it
					 */
						blist_append(layer[n].b2gmap[i][j], u, 1.0);
					//	warning("overlap of functional blocks?\n");
					}
				/* boundary grid cells partially overlapped by this unit	*/
				} else {
					/* shortcuts for cell boundaries	*/
					double lc = j * cw, rc = (j+1) * cw;
					double tc = stk->length - i * ch;
					double bc = stk->length - (i+1) * ch;

					/* shortcuts for overlap width and height	*/
					double oh = (MIN(tu, tc) - MAX(bu, bc));
					double ow = (MIN(ru, rc) - MAX(lu, lc));
					double occupancy;

					/* overlap tolerance	*/
					if (eq(oh/ch, 0))
						oh = 0;
					else if (eq(oh/ch, 1))
						oh = ch;

					if (eq(ow/cw, 0))
						ow = 0;
					else if (eq(ow/cw, 1))
						ow = cw;

					occupancy = (oh * ow) / (ch * cw);
					if (oh < 0 || ow < 0)
						fatal("negative overlap!\n");

					/* first unit in the list	*/
					if (!layer[n].b2gmap[i][j])
						layer[n].b2gmap[i][j] = new_blist(u, occupancy);
					else
					/* append at the end	*/
						blist_append(layer[n].b2gmap[i][j], u, occupancy);
				}
	}

	/*
	 * sanity check
	test_b2gmap(model, layer);
	 */
}


/*calculate total number of units of the 3D stack*/
int obtain_total_units(stk_descr_t *stk)
{
	int n_total_units=0;
	int i;

	for(i=0;i<stk->n_layers;i++)
	{
		n_total_units+=stk->layers[i].flp->n_units;
	}

	return n_total_units;
}

/*debug use, print out stack desciption*/
void check_stk_description(stk_descr_t* stk)
{
	printf("chip width: %f\n", stk->width);
	printf("chip length: %f\n", stk->length);
	printf("number of die layers: %d\n", stk->n_layers);
	printf("grid rows: %d, cols: %d\n", stk->rows, stk->cols);

	if(stk->cooling_method==1)
	{
		printf("heatsink thickness: %f; width(length): %f; conductivity: %f; R_conv: %f\n",
				stk->heatsink->t_heatsink, stk->heatsink->width_heatsink, stk->heatsink->k_heatsink, stk->heatsink->r_conv);
		printf("spreader thickness: %f; width(length): %f; conductivity: %f\n",
				stk->heatsink->t_spreader, stk->heatsink->width_spreader, stk->heatsink->k_spreader);
		printf("TIM thickness: %f; conductivity: %f\n", stk->heatsink->t_TIM, stk->heatsink->k_TIM);
	}
	else if (stk->cooling_method==2)
	{
		printf("channel width: %f\n", stk->microchannel->channel_width);
		printf("channel height: %f\n", stk->microchannel->channel_height);
		printf("wall width: %f\n", stk->microchannel->wall_width);
		printf("number of channels per layer: %d\n", stk->microchannel->n_channels_per_layer);
		printf("default total flow rate: %f\n", stk->microchannel->default_flowrate);
	}

	int i;

	for(i=0;i<stk->n_layers;i++)
	{
		printf("layer %d:\tBEOL %e\tSource %e\tBulk %e\n", stk->layers[i].idx, stk->layers[i].BEOL_height,stk->layers[i].source_height,stk->layers[i].bulk_height);
	}

}

void tsv_translate(stk_descr_t* stk)
{
	double cw = stk->width/stk->cols;
    double ch = stk->length/stk->rows;
    int nr=stk->rows;
    int nc=stk->cols;
    tsv_physi_t *phy=stk->tsv_physi;
    int i;
    int n_ch_per_cell=stk->microchannel->n_channels_per_layer/stk->cols;


    init_tsv(stk);

    double west_tsv, east_tsv, north_tsv, south_tsv;
    double lc, rc, tc, bc, oh, ow;
    double kxy_tsv, kz_tsv, dm;
    double gamma=stk->k_tsv/K_Liner, betaxy, betaz;
    int xleft, xright, ybottom, ytop;
    int n, j, k, layer, sublayer;
    double xcenter, ww;
    int wallindex;

    for(n=0;n<stk->n_tsvs;n++)
    {
    	dm=phy[n].width;
    	kxy_tsv=(1+dm*dm*(gamma-1)/((dm+2*T_Liner)*(dm+2*T_Liner*gamma)))*K_Liner;
    	kz_tsv=(1+dm*dm*(gamma-1)/((dm+2*T_Liner)*(dm+2*T_Liner)))*K_Liner;

    	west_tsv=phy[n].x-0.5*phy[n].width-T_Liner;
    	east_tsv=phy[n].x+0.5*phy[n].width+T_Liner;
    	north_tsv=phy[n].y+0.5*phy[n].width+T_Liner;
    	south_tsv=phy[n].y-0.5*phy[n].width-T_Liner;

    	xleft=double2int(west_tsv/cw);
    	xright=double2int(east_tsv/cw);
    	ybottom=nr-double2int(north_tsv/ch)-1;
    	ytop=nr-double2int(south_tsv/ch)-1;

    	for(i=xleft;i<xright+1;i++)
    		for(j=ybottom;j<ytop+1;j++)
    		{
				/* shortcuts for cell boundaries	*/
    			lc = i * cw;
    			rc = (i+1) * cw;
				tc = stk->length - j * ch;
				bc = stk->length - (j+1) * ch;

				/* shortcuts for overlap width and height	*/
				oh = (MIN(north_tsv, tc) - MAX(south_tsv, bc));
				ow = (MIN(east_tsv, rc) - MAX(west_tsv, lc));

				/* overlap tolerance	*/
				if (eq(oh/ch, 0))
					oh = 0;
				else if (eq(oh/ch, 1))
					oh = ch;

				if (eq(ow/cw, 0))
					ow = 0;
				else if (eq(ow/cw, 1))
					ow = cw;

				if (oh < 0 || ow < 0)
					fatal("negative overlap of TSV and grid cells!\n");

				for(sublayer=phy->z1*NumSubLayers+phy->sl1;sublayer<phy->z2*NumSubLayers+phy->sl2+1;sublayer++)
				{
					if(sublayer%3==0)
					{
						betaxy=kxy_tsv/K_BEOL;
						betaz=kz_tsv/K_BEOL;
					}
					else
					{
						betaxy=kxy_tsv/K_SILICON;
						betaz=kz_tsv/K_SILICON;
					}
					stk->tsv_circuitx[nr*nc*sublayer+j*nc+i]+=ow*oh*(betaxy-1)/(cw*ch*(1+(1-ow/cw)*(betaxy-1)));
					stk->tsv_circuity[nr*nc*sublayer+j*nc+i]+=ow*oh*(betaxy-1)/(cw*ch*(1+(1-oh/ch)*(betaxy-1)));
					stk->tsv_circuitz[nr*nc*sublayer+j*nc+i]+=ow*oh*(betaz-1)/(cw*ch);
				}


				if(stk->cooling_method==1)
				{
					betaxy=kxy_tsv/stk->heatsink->k_TIM;
					betaz=kz_tsv/stk->heatsink->k_TIM;

					for(sublayer=phy->z1;sublayer<phy->z2;sublayer++)
					{
						stk->tsv_bondx[nr*nc*sublayer+j*nc+i]+=ow*oh*(betaxy-1)/(cw*ch*(1+(1-ow/cw)*(betaxy-1)));
						stk->tsv_bondy[nr*nc*sublayer+j*nc+i]+=ow*oh*(betaxy-1)/(cw*ch*(1+(1-oh/ch)*(betaxy-1)));
						stk->tsv_bondz[nr*nc*sublayer+j*nc+i]+=ow*oh*(betaz-1)/(cw*ch);
					}

				}
				else if(stk->cooling_method==2)
				{
					xcenter=0.5*(MIN(east_tsv, rc) + MAX(west_tsv, lc))-lc;
					wallindex=double2int((double2int(2*xcenter/(stk->microchannel->channel_width+stk->microchannel->wall_width))+1)/2);
					betaxy=kxy_tsv/K_SILICON;
					betaz=kz_tsv/K_SILICON;
					if(wallindex==0||wallindex==n_ch_per_cell)
						ww=0.5*stk->microchannel->wall_width;
					else
						ww=stk->microchannel->wall_width;
					for(sublayer=phy->z1;sublayer<phy->z2;sublayer++)
					{
						stk->tsv_bondx[(nr*nc*sublayer+j*nc+i)*(n_ch_per_cell+1)+wallindex]+=ow*oh*(betaxy-1)/(ww*ch*(1+(1-ow/ww)*(betaxy-1)));
						stk->tsv_bondy[(nr*nc*sublayer+j*nc+i)*(n_ch_per_cell+1)+wallindex]+=ow*oh*(betaxy-1)/(ww*ch*(1+(1-oh/ch)*(betaxy-1)));
						stk->tsv_bondz[nr*nc*sublayer+j*nc+i]+=ow*oh*(betaz-1)/(ww*ch);
					}
				}
    		}
    }
}

void init_tsv(stk_descr_t* stk)
{
	int nl=stk->n_layers;
    int nr=stk->rows;
    int nc=stk->cols;
    int i;


    /*allocate space for TSV map on circuit layers in grid level*/
    stk->tsv_circuitx = (double*)calloc(nr*nc*nl*NumSubLayers, sizeof(double));
    stk->tsv_circuity = (double*)calloc(nr*nc*nl*NumSubLayers, sizeof(double));
    stk->tsv_circuitz = (double*)calloc(nr*nc*nl*NumSubLayers, sizeof(double));

    for(i=0;i<nr*nc*nl*NumSubLayers;i++)
    {
    	 stk->tsv_circuitx[i]=0;
    	 stk->tsv_circuity[i]=0;
    	 stk->tsv_circuitz[i]=0;
    }

    /*allocate space for TSV map on bonding layers, i.e. TIM layers (heatsink cooling), or channel layers (microfluidic cooling)*/
    if(stk->cooling_method==1)
    {
    	stk->tsv_bondx=(double*)calloc(nr*nc*(nl-1), sizeof(double));
    	stk->tsv_bondy=(double*)calloc(nr*nc*(nl-1), sizeof(double));
    	stk->tsv_bondz=(double*)calloc(nr*nc*(nl-1), sizeof(double));
    	for(i=0;i<nr*nc*(nl-1);i++)
    	{
    		stk->tsv_bondx[i]=0;
    		stk->tsv_bondy[i]=0;
    		stk->tsv_bondz[i]=0;
    	}

    }
    else if(stk->cooling_method==2)
    {
    	int n_ch_per_cell=stk->microchannel->n_channels_per_layer/stk->cols;
    	stk->tsv_bondx=(double*)calloc(nr*nc*(nl-1)*(n_ch_per_cell+1), sizeof(double));
    	stk->tsv_bondy=(double*)calloc(nr*nc*(nl-1)*(n_ch_per_cell+1), sizeof(double));
    	stk->tsv_bondz=(double*)calloc(nr*nc*(nl-1)*(n_ch_per_cell+1), sizeof(double));
    	for(i=0;i<nr*nc* (nl-1)*(n_ch_per_cell+1);i++)
    	{
    		stk->tsv_bondx[i]=0;
    		stk->tsv_bondy[i]=0;
    		stk->tsv_bondz[i]=0;
    	}
    }
}
